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SUMMARY 


This  report  describee  progress  made  in  research  during  the  final  year  of 
a  five-year  study  of  numerical  dynamical  analysis  of  structures  funded  by  an 
AFOSR  Grant.  The  original  proposal  was  JPL  Proposal  51-641  dated  6-29-76, 
and  the  final  year  renewal  was  proposed  on  12-18-80. 

The  proposed  work  effort  in  the  final  year's  work  Included  the  following 
Tasks: 

Task  1  -  Development  of  a  higher  order  rectangular  plane  stress/strain 

finite  element. 

Task  II  -  Development  of  a  solid  hexahedron  finite  dynamic  element. 

Task  III  -  Further  refinement  of  the  Associated  Generalized  Eigenproblem 
Solution  Routine. 

The  nature  of  this  work  is  to  generate  analytical  results  for  publication 
in  the  open  literature.  It  is  connton  procedure,  therefore,  to  report  the 
results  of  contractural  work  by  submitting  preprints  or  reprints  of  articles 
to  be  published  as  the  result  of  research  supported  under  these  tasks.  We 
have  therefore  collected  together  the  appropriate  preprints  and  reprints  and 
packaged  them  together  as  the  report  for  the  final  year's  effort. 

In  addition  to  this  work  in  strict  adherance  with  the  Task  descriptions, 
a  piece  of  research  was  carried  out  on  the  application  of  the  finite  dynamic 
element  method  to  coupled  fluid-structure  problems.  This  was  done  at  the 
Invitation  of  the  Fourth  International  Symposium  on  Finite  Element  Methods  in 
Flow  Problems,  and  represents  a  logical  extension  of  developments  within  the 
AFOSR  Grant  to  new  application  areas  of  interest  to  the  Air  Force. 

The  report  is  therefore  divided  into  four  parts: 

Task  I  -  Higher  order  element 

Task  II  -  Solid  hexahedron  element 

Task  III  -  Refinement  of  the  solution  routine 

Task  IV  -  Application  tp  fluid-structure  problems 


REPORT  OR  TASK  I 


"Development  of  a  Higher-Order  Plane  Finite 
Dynamic  Element" 


(to  be  published) 


NUMERICAL  FORMULATION  FOR  A  HIGHER-ORDER 
PLANE  FINITE  DYNAMIC  ELEMENT 


K.  K.  Gupta 


A. 
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SUMMARY 


The  paper  is  essentially  concerned  with  the  development  of  a  8-node  plane 
rectangular  finite  dynamic  element  and  presents  detailed  descriptions  of  the 
associated  numerical  formulation  involving  the  higher  order  dynamic  correction 
terms,  pertainings  to  the  related  stiffness  and  inertia  matrices. 

Numerical  test  results  of  free  vibration  analysis  are  presented  in  detail 
for  the  newly  developed  8-node  element,  as  well  as  the  corresponding  4-node 
element  to  afford  a  clear  comparison  of  the  relative  efficiences  of  the  corres* 
ponding  finite  element  and  the  dynamic  element  procedures.  Such  results  indi¬ 
cate  a  superior  pattern  of  solution  convergence  of  the  presently  developed 
dynamic  element. 


INTRODUCTION 


A  discrete  element  Idealization  of  a  continuum,  undergoing  free  vibration, 
may  be  achieved  by  uniquely  relating  the  displacement  field  within  an  element  in 
terms  of  its  nodal  values.  Such  a  relationship  is  expressed  as 

jj  ■  a(w)'  IJ  (1) 

in  which  the  shape  function  matrix  a  is  a  function  of  the  natural  frequencies  u, 
of  the  structure  under  consideration.  The  associated  stiffness  and  Inertia 
matrices  may  then  be  derived  by  standard  procedures,  based  on  variational  princi¬ 
ples,  noting  that  the  resulting  matrices  are  obtained  as  functions  of  the  initially 
unknown  natural  frequencies.  Subsequent  extraction  of  roots  and  vectors  from 
these  matrices  is  extremely  difficult  and  uneconomical  in  nature  and  to  avoid 
such  a  unwieldy  formulation  equation  (1)  Is  expanded  in  ascending  powers  of  u>, 
as  below: 
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ji(w)  ■  Sq  +  ua^  +  “  *2  +  *•*  (2) 

resulting  In  the  following  expressions  for  the  element  stiffness  and  Inertia 
matrices 

K  -  Kq  +  +  w4^  +  . . .  (3) 

M  -  Mq  +  +  o)4^  +  . . .  (4) 

These  matrices  when  appropriately  combined  yields  the  global  matrices  for  the 
entire  continuum.  The  associated  free  vibration  formulation  is  given  by 

[Kq-^CMq-^)  -  -...]£«  0  (5) 


in  which  the  higher  order  dynamic  correction  matrices  K^,  K^,  Mj  and  such  other 
terms  are  retained  in  the  formulation  in  the  dynamic  element  method  (DEM) ,  whereas 
only  the  Initial  terms  Kq  and  Mq  are  included  in  the  analysis  employing  the  usual 
finite  element  method  (FEM).  Furthermore,  in  the  dynamic  element  procedure,  the 
series  form  of  equation  (3)  is  suitably  truncated  to  yield  a  quadratic  matrix 
eigenvalue  problem 


(A-XB-X  C)  _g  -  0  (6) 

2 

with  A  -  Kq,  £  -  Mq-K^ ,  £  -  M^-K^  and  X  -  to  ,  whereas  in  the  finite  element  method 
the  equivalent  formulation  has  the  form 

(A-XB)  -  0  (7) 


where 


1“  Mq  * 


A  dynamic  element  formulation  was  earlier  achieved  for.  line  elements  ,  involv 
ing  derivation  of  the  higher  order  dynamic  correction  terms.  The  procedure  was 
further  developed  for  continuum  discretization  and  details  of  the  relevant  dynamic 
elements  pertaining  to  membrane  and  plane  stress  *  problems  have  been  published 
earlier.  Such  results  Indicate  that  dynamic  elements  exhibit  much  superior  solu¬ 
tion  convergence  characteristics  when  compared  with  the  usual  finite  element 
method,  resulting  in  substantial  economy  in  the  free  vibration  analysis  of  prac¬ 
tical  problems.  Furthermore,  the  usual  solution  techniques5  for  the  quadratic 
matrix  equation  Involved  the  eigenproblem  solution  of  an  equivalent  system 
characterized  by  a  single  full  matrix  of  order  twice  that  of  the  original  system, 
requiring  prohibitive  computational  effort  forjaost  practical  problems.  However, 
new  solution  techniques***^  for  the  quadratic  matrix  formulation  of  equation  (6) 
pertaining  to  the  dynamic  element  method  enable  eigenproblem  solution  with 
approximately  the  same  computational  effort  as  that  required  for  the  solution  of 
equation  (7)  associated  with  a  finite  element  formulation. 

The  purpose  of  this  paper  is  to  present  detailed  formulation  of  a  8-node, 
plane  rectangular  dynamic  element,  numerical  results  are  presented  for  a  repre¬ 
sentative  problem,  solved  by  both  the  DEM  and  FEM  formulation.  Furthermore, 
similar  results  are  presented  for  the  corresponding  4-node  element  to  afford  a 
clear  comparison  of  convergence  characteristics  of  the  various  element  types. 


DYNAMIC  ELEMENT  FORMULATION  FOR  A  PLANE 
8-NODE  RECTANGULAR  ELEMENT 


Figure  1  depicts  a  typical  rectangular,  8-node  plane  element.  The  differen¬ 
tial  equations  of  free  vibration  of  such  a  continuum  are  of  the  following  form 
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where  p,  v>  and  E  are  the  element  mass  per  unit  area,  Poisson's  ratio  and  the 
Young's  modulus,  respectively.  Appropriate  solution  of  equations  (8)  and  (9) 
for  the  in-plane  deformations  ux,  u^  may  be  expressed  ir  infinite  series  form 
as  below 


u  ■  a(u)U  »  (a-  +  ua,  +  u>  a_  +  ...)  U 

y  —  —  —Ox  —lx  — 2x  — 


uy  *  «<“>£  *  ^oy  +  +  “  a2y  +  '  *  0  - 


(10) 
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which  if  substituted  in  equations  (8)  and  (9)  yields  the  final  expressions  for 
the  differential  equations.  As  for  example,  such  equations  in  the  y-direction  are 
as  follows 


.2  ,2  .2  J2. 
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in  which 


“l  ’  (l-2y)  ’  a2  *  (l-2y)  and  6  =  V  (1+u) ’ 


similar  expressions  being  also  obtained  in  the  x  direction. 

When  solving  such  differential  equations,  it  is  postulated  that  while  a^  is 
allowed  to  satisfy  the  appropriate  boundary  conditions  a^  and  must  all  vanish 
at  the  boundaries.  Thus  the  solutions  for  a^x*  a^y  satisfying  equation  (1)  and 
its  counterpart  are  assumed  as  below 

^Ox  *  C1  +  C2X  +  C3y  +  C4x2  +  C5Xy  +  c6y2  +  c7x2y  +  C8Xy2  (15) 

-^Oy  c9  +  C10X  +  Clly  +  C12X  +  C13Xy  +  C14y  +  C15X  y  +  C16Xy 


in  which  the  coefficients  c^  -  c^g  are  evaluated  by  satisfaction  of  the  boundary 
conditions 


ux  -  Ux,  uy  -  U2  at  x  -  0,  y  -  0  u^  -  Uy  uy  -  at  x  ■  2a,  y  -  0 
\  -  U5  uy  -  Ug  at  x  -  2a,  y  -  2b  uv  -  U?,  -  UR  at  x  -  0,  y  -  2b 


x  '  y  8 


Ux  "  U9  uy  '  U10  at  x  •  a,  y  -  0  Ujc  -  Uu,  uy  -  U12  at  x  -  2a,  y  -  b 

Ux  “  U13  Uy  "  Ui4  at  x  "  a»  7  "  2b  ux  "  U15»  uy  "  ul6  at  x  ■  °»  y  *  b 

Expressions,  similar  to  equations  (15)  and  (16)  are  chosen  as  solutions  for  a., 
yielding  a^x  ■  J),  a^y  ■  j);  appropriate  solutions  for  equation  (14)  and  its  counter 
part,  on  the  other  hand,  are  assumed  as  follows 


~2x  ■  S1  +  d2*  +  %  +  V*  +  V*  +  V  +  +  d3Xy2 
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when  the  coefficients  of  the  complementary  functions  are  determined  by 

satisfying  the  appropriate  boundary  conditions,  noting  that  the  coefficients 
c^  -  c^g  are  computed  earlier  from  equations  (15)  and  (16).  The  final  forms  of 
the  shape  functions  are  obtained  after  performing  some  routine  algebraic  manipu¬ 
lations.  Appropriate  derivations  of  the  higher  order  shape  function  matrix  a^ 
is  crucial  in  the  development  of  the  current  finite  dynamic  element. 

The  shape  function  vectors  are  thus  defined  as 

a  ■  a,,  +  u>2a„  »  a  -  a..  +  u)2a-  (19) 

—x  -^Ox  — 2x  — y  — Oy  — 2y 

in  which  the  scalars  of  each  vector  are  coupled  to  the  appropriate  nodal  degrees 
of  freedom  of  the  element.  Associated  strain-displacement  relationship  is  given 

as 


1-8 


(20) 


in  which 
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The  individual  strain-displacement  matrices  are  suitably  combined  to  yield 
the  bQ  and  matrices. 

The  element  stiffness  and  inertia  matrices  may  next  be  developed  by  standard 
procedure  once  the  various  a  and  t>  matrices  have  been  determined,  as  above.  Thus 
the  stiffness  matrix  is  obtained  as 


K  -  Kq  +  u) 


where 
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»  V 

■/  V*  it 
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and  in  which  x  is  the  stress-strain  matrix  for  two  dimensional  elasticity,  v 
being  volume  of  the  element.  In  a  similar  manner  the  inertia  matrices  are  given 


s*  -  So,  +  ■  Sj* 


m  -  +  oi 


S0y  +  «  S2y 


.V-V* 


in  which 


*0.  ■  °v/ ioH  4o*4v 

*2*  '  °,f  4o»I*2*dv 

S'’./Vs or4" 

■2,  *  »u/loT%y4’ 

where  0  La  tha  mesa  par  wait  volume  of  cha  plata  alaaant. 

The  symbolic  manipulation  program  MACSYMA® ,  has  baan  utilized  in  processing 
aquations  (15)-(30)  for  tha  derivation  of  tha  a,  t>,  m  and  K  matrices  Involving 
rathar  large  amaamta  of  algebraic  aaaipulatlons.  The  resulting  expression  for 
tha  matrix  elamaata  era  next  traaaformad  in  FORTRAN  prograamed  form,  by  employing 
a  suitable  MhCSTMA  laat ruction.  Dua  to  tha  lengthy  nature  of  the  expressions  for 
£»  •»  K  matrices,  they  are  aot  reproduced  hare  for  ready  reference; 

however,  tha  programmed  form  of  these  individual  expressions  may  be  supplied  for 
appropriate  utilisation.  Tha  element  matrices  are  than  combined  ty  standard 
process  to  yield  tha  global  stiffness,  inertia  and  dynamic  correction  matrices 
for  subsequent  analysis  of  numerical  examples,  by  solving  the  quadratic  eigenvalue 
problem  depicted  by  equation  6. 


(27) 
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v  I  — 2x  —Ox 
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NUMERICAL  RESULTS 


A  square  plate,  (Figure  2)  with  one  edge  fixed  and  three  edges  free,  vibrating 

3  4 

in  its  own  plane,  is  chosen  as  a  numerical  example  as  before  '  ,  having  the 
following  basic  structural  data:  side  length  (1)  -  10.0,  thickness  (t)  -  1.0; 
Poisson's  ratio  (y)  -  0.3.  Solution  results  were  obtained  for  the  model  employ¬ 
ing  an  increasing  number  of  elements  and  such  analyses  were  performed  for  both 
the  dynamic  and  the  finite  element  idealizations.  Table  I  presents  these  results 
in  parametric  form,  along  with  similar  results  obtained  by  utilizing  4-node 

3 

rectangular  elements  .  Such  a  table,  on  the  other  hand,  provides  a  clear  com¬ 
parison  of  the  pattern  of  root  convergence  of  the  higher-order  8-node  and  the 
simple  4-node  dynamic  and  finite  elements,  in  a  concised  form. 

Figure  3  depicts  the  pattern  of  convergence  of  two  typical  roots  pertaining 
to  the  four  sets  of  results,  for  varying  mesh  sizes.  Such  results  are  also 
depicted  in  Figure  4,  as  a  function  of  total  computational  effort  involved  in  the 
respective  eigenproblem  solution.  The  solution  results  pertaining  to  a  20  x  20 
mesh  discretization  is  accepted  as  the  exact  solution,  in  the  absence  of  an 
available  analytical  solution. 


CONCLUDING  REMARKS 


Numerical  results  presented  in  Table  I  are  further  depicted  in  Figures  3 
and  4  to  provide  a  better  insight  into  the  patterns  of  behavior  of  the  newly 
developed  8-node  dynamic  element  as  veil  as  its  finite  element  counterpart.  Such 
results  are  also  presented  for  the  simple  4-node  element  to  afford  a  clear  com¬ 
parison  of  the  relative  solution  efficiencies  for  both  the  elements  employing 
either  the  dynamic  element  or  the  finite  element  technique.  It  is  quite  apparent 
from  these  results  that  significant  improvement  in  root  convergence  is  achieved 
when  dynamic  elements  are  used  in  place  of  the  usual  finite  elements.  Furthermore, 
Figure  4  indicates  that  a  4-node  dynamic  element  displays  convergence  characteris¬ 
tics  similar  to  an  8-node  finite  element.  Thus,  for  a  required  two  percent 
solution  accuracy  for  the  elgenproblem  solution  efforts  for  the  8-node  DEM/FEM 
and  4-node  DEM/FEM  procedures  bear  the  ratios  1,6,4  and  15  respectively.  Also, 
with  Increasing  mesh  size,  errors  in  frequencies  computed  by  the  DEM  analysis 
decrease  much  more  rapidly  than  the  FEM  computations.  Furthermore,  for  a  given 
solution  accuracy,  the  DEM  analysis  requires  considerably  less  data  preparation 
effort  due  to  a  significant  reduction  in  mesh  size. 

As  pointed  out  in  the  Introduction  section,  the  development  of  the  dynamic 

elements  proved  to  be  highly  beneficial  only  after  new  elgenproblem  solution 

techniques  were  formulated  that  enable  solution  of  the  quadratic  matrix  equation 
2  2 

(A  “  <*>  A  ~  u  G)  jl  ■  0.  A  discussion  on  the  choice  of  the  higher  order  shape 

4- 

function,  which  is  crucial  to  the  current  formulation,  is  given  elsewhere  . 
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Table  I 


Comparison  of  natural  frequency  values  of  a  square  cantilever  plate  (FEM  and  OEM 
results  for  8-node  elements  shown  in  upper  and  lower  rows,  corresponding  values 
for  4-node  elements  shown  in  parenthesis) 


Eigenvalue  Parameter  u  ■ 

ui >/  (E/P) 

Mesh 

A 

Size 

"l 

“2 

“3 

“4 

“5 

“6 

lxl 

.07046 

(.07792) 

.1602 

(.1743) 

.1926 

(.2908) 

.3169 

.3867 

.3968 

.07027 

(.07444) 

.1522 

(.1491) 

.1736 

(.2444) 

.2715 

.3182 

.3708 

2x2 

.06708 

.1585 

.1817 

.2906 

.3167 

.3289 

(.07186) 

(.1637) 

(.2090) 

(.3372) 

(.3905) 

(.3964) 

.06706 

.1579 

.1797 

.2801 

.3058 

.3077 

(.07096) 

(.1547) 

(.1946) 

(.2960) 

(.3340 

(.3441) 

3x3 

.0638 

.1583 

.1785 

.2835 

.3083 

.3235 

(.06913) 

(.1608) 

(.1934) 

(.3152) 

(.3500) 

(.3609) 

.06636 

.1582 

.1781 

.2803 

.3049 

.3181 

(.06876) 

(.1565) 

(.1867) 

(.2923) 

(.3190) 

(.3288) 

20x20 

(exact 

results) 

(.06585) 

(.1579) 

(.1769) 

(.2796) 

(.3033) 

(.3214) 
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ERROR  IN  FREQUENCY, 


Figure  4.  Plane  rectangular  elements:  Comparison  of  total  numerical  efforts  as 
a  function  of  error  in  frequencies 
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REPORT  ON  TASK  II 


"Development  of  a  Solid  Hexahedron 
Finite  Dynamic  Element" 


(to  be  published  as  a  full  paper) 


DEVELOPMENT  OF  A  SOLID  HEXAHEDRON  FINITE  DYNAMIC  ELEMENT 

A  solid  hexahedron  element,  with  three  translational  degrees  of  freedom  per 
node,  is  shown  in  Figure  1.  The  differential  equations  of  free  vibrations  of  the 
associated  continuum  may  be  expressed  as 
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U,p,E  being  the  parameters  as  defined,  earlier.  Solutions  of  Eqs.  (1),  (2)  and  (3) 
are  taken  respectively,  as 
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(10) 


y  —  Direction 
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Solutions  of  Equations  (8) ,  (11)  and  (14)  are  assumed  in  the  following  from 


Ux  -  aQjcU  ■  +  C2»  +  Cgy  +  C^z  +  C5xy  +  Cgyz  +  C^xz  +  Cgxyz  (17 

uy  "  a0y°  "  C9  +  Ciox  +  CU?  +  C12Z  +  Cl3xy  +  Cl^z  +  ci5xz  +  <18 

Uz  *  a0zU  “  C17  +  C18x  +  C19y  +  C20z  +  Sl3^  +  C22yz  +  C23XZ  +  C2^z  (19 


k* 


where  the  coefficients  -  C ^  are  evaluated  by  appropriate  satisfaction  of  the 
boundary  conditions  as  below: 
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Once  the  shape  function  matrices  and  the  strain-displacement  relations  have  been 
derived  the  mass  and  stiffness  matrices  may  be  formulated  in  a  manner,  similar  to 
the  preceding  section.  As  before,  the  MACSYMA  program  has  been  extensively  used 
for  all  algebraic  manipulations  involved  in  the  development  of  the  various  a,  b, 
m  and  fc  matrices. 


NUMERICAL  EXAMPLE 


A  cub*  of  unit  dimensions  i*  chosen  as  a  numerical  example.  Both  finite 
and  dynamic  element  discretizations  were  achieved  by  varying  mesh  sizes. 
Results  of  such  analyses  are  given  in  Table  I. 

Table  I.  Natural  frequencies  of  a  cube  (FEM  and  DEM  results  are 
shown  in  upper  and  lower  rows  respectively) 


Mesh 

Size 


Eigenvalue  Parameters 


.7928 

.7361 

.7930 

.7361 

1.0742 

.9883 

1.7845 

1.3993 

2.9197 

2.1052 

2.9417 

2.1052 

.7378 

.7228 

.7381 

.7228 

.9994 

.9666 

1.6947 

1.5492 

2.1001 

1.8761 

2.1003 

1.8762 

.7077 

.7015 

.7078 

.7015 

.9598 

.9411 

1.6478 

1.5777 

1.9388 

1.8348 

1.9389 

1.8348 

The  above  results  indicate  that  the  dynamic  elements  are  significantly  more 
efficient  than  the  corresponding  finite  elements. 
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report  on  task  III 


"Development  of  a  Unified  Numerical  Procedure  for  Free  Vibration 
Analysis  of  Structures",  IJNME,  1981. 


"STARS  —  Structural  Analysis  Routines",  Data  Input  Procedure. 


STARS  (Structural  Analysis  Routines) 
K.K. Gupta 


DATA  INPUT  PROCEDURE 


PROGRAM  DESCRIPTION 


Stars  is  a  compact,  general-purpose,  finite  and  dynamic 
element  computer  program  for  static  (multiple  load) ,  free 
vibration  and  dynamic  response  analysis  of  undamped  and 
damped  (viscous  and  structural)  structures  including 
spinning  ones. 

The  element  library  presently  consists  of  line  and 
triangular/quadrilateral  shell  elements  enabling  solution  of 
truss,  beam,  space  frame,  plane,  plate,  shell  structures,  or 
any  combination  thereof.  The  program  allows  zero,  finite  and 
interdependent  deflection  boundary  conditions.  The  dynamic 
response  analysis  capability  allows  initial  deformation  and 
velocity  inputs  whereas  the  transient  excitation  may  be  either 
force  of  acceleration  in  nature. 

Data  input  can  be  at  random  within  a  specific  data  set, 
and  the  program  offers  certain  automatic  data  generation 
capabilities.  Input  data  is  formated  as  an  optimal  combination 
of  free  and  fixed  formats. 

The  program,  developed  in  modular  form  for  easy  modifications 
is  written  in  FORTRAN  for  the  VAX  11/780  computer  and  consists  of 
about  6000  instructions.  Continued  development  of  the  program  is 
envisaged  while  keeping  tight  control  on  the  size  of  the  program. 
Extensive  interactive  plotting  capabilities  form  an  important 
feature  of  the  program. 


DESCRIPTION  OF  MAIN IB  INPUT 


CARD 


CARD 


1  FORMAT 
JOB  TITLE 
13A6 


2  FORMAT 

NN,NfNEL,NMAT,NMECN,NEP,NET,NTMP,NPR,NBUN 

FREE  FORMAT 

WHERE: 

NN  -  NUMBER  OF  NODES 
N  -  ORDER  OF  PROBLEM  (USUALLY  6*NN) 

NEL  -  TOTAL  NUMBER  OF  ELEMENTS 

NMAT  -  TOTAL  NUMBER  OF  ELEMENT  MATERIAL  TYPES 

NMECN  -  NUMBER  OF  MATERIAL  ELASTIC  CONSTANTS 

NEP  -  TOTAL  NUMBER  OF  ELEMENT  PROPERTY  TYPES  (LINE  ELEMENTS) 
NET  -  TOTAL  NUMBER  OF  ELEMENT  THICKNESS  TYPES  (SHELL  ELEMENTS) 
NTMP  -  TOTAL  NUMBER  OF  ELEMENT  TEMPERATURE  TYPES 
NPR  -  TOTAL  NUMBER  OF  ELEMENT  PRESSURE  TYPES 
NBUN  -  NUMBER  OF  NODAL  CONNECTIVITY  CONDITIONS 


3-3 


CARD  3  FORMAT 

1PROB , I BAN , NPREC ,NC , I DRS ,  IPLOT 

FREE  FORMAT 

WHERE: 

IPROB  -  INDEX  FOR  PROBLEM  TYPES,  TO  BE  SET  AS  FOLLOWS: 

-  1,  UNDAMPED  FREE  VIBRATION  CASE  FOR  NON -SPINNING 

STRUCTURES. 

-  2,  UNDAMPED  FREE  VIBRATION  CASE  FOR  SPINNING  STRUCTURES, 

C  BEING  SKEW -SYMMETRIC  CORIOLIS  MATRIX,  K-KE+KG+KC 
WHERE  KE , KG , KC  ARE  RESPECTIVELY  THE  ELASTIC, 
GEOMETRICAL  STIFFNESS  AND  CENTRIFUGAL  FORCES  MATRIX 
RESPECTIVELY 

-  3,  QUADRATIC  MATRIX  El GEN PROBLEM  OPTION  FOR  DEM  ANALYSIS, 

K-LAMBDA*M-LAMBDA**2*C,  C  IS  DYNAMIC  CORRECTION 
MATRIX. 

-4,  FOR  SPINNING  STRUCTURES,  WHEN  C-CC+CD,  CC  IS  THE 
SKEW -SYMMETRIC  CORIOLIS  MATRIX,  CD  BEING  THE 
DIAGONAL  VISCOUS  DAMPING  MATRIX  IN  R-I*LAMBDA*C 
-LAMBDA**2*M  FORMULATION. 

-  5,  AS  IN  4,  WITH  STRUCTURAL  DAMPING,  K  BEING  COMPLEX 

-  6,  FOR  NON-SPINNING  STRUCTURES  WITH  VISCOUS 

DAMPING  (C). 

-  7,  AS  IN  6,  WITH  STRUCTURAL  DAMPING 

-  8,  SOLVES  SIMULTANEOUS  EQUATIONS  AX-B  ONLY,  WHEN 

MATRIX  A,  REAL  SYMMETRIC  OR  HERMITIAN,  IS  READ  IN 
LIEU  OF  K. 

IBAN  -  BANDWIDTH  MINIMIZATION  OPTION 
»  0,  PERFORM  MINIMIZATION 

-  1,  DO  NOT  PERFORM  MINIMIZATION 
NPREC  -  PRECISION  FOR  SOLVING  PROBLEM 

-  1,  SINGLE  PRECISION  REAL 

-  2.  DOUBLE  PRECISION  REAL 

-  3,  SINGLE  PRECISION  COMPLEX  * 

-  4,  DOUBLE  PRECISION  COMPLEX 

NC  -  NUMBER  OF  SETS  OF  CONCENTRATED  NODAL  LOAD/MASS  DATA 
(IF  NONE  SET  TO  1) 

IDRS  -  INDEX  FOR  DYNAMIC  RESPONSE  ANALYSIS 

-  0,  NO  RESPONSE  ANALYSIS  NEEDED 

-  1,  PERFORM  RESPONSE  ANALYSIS 
IPLOT  -  INDEX  FOR  PLOTTING 

-  0,  NO  PLOTS  ARE  REQUIRED 

-  1,  PLOT  STRUCTURE  GEOMETRY;  RESTART  TO  CONTINUE  SOLUTION 

-  2,  PLOT  FINAL  NODAL  DEFORMATIONS  AS  FUNCTIONS  OF  TIME, 
MODE  SHAPES  AND  ELEMENT  STRESSES 


CARD  4  FORMAT 

IMLUMP , I PLUMP , I SPIN , I PRINT 
FREE  FORMAT 
WHERE • 

IMLUMP  -  INDEX  FOR  NODAL  LUMPED  MASS  INPUT 

-  0  ,  NO  LUMPED  MASS  INPUT 

-  1  ,  LUMPED  MASS  DATA  INPUT  (IPROB-1  THRU  7) 

I PLUMP  -  INDEX  FOR  NODAL  EXTERNAL  LOAD  INPUT 

-  0  ,  NO  LOAD  INPUT 

-  1  ,  CONCENTRATED  NODAL  LOAD  INPUT  (IPROB-8) 

ISPIN  -  SPIN  OPTION  FOR  STRUCTURE 

-  0,  NO  SPIN 

-  1 ,  SPINNING  STRUCTURE 
I PRINT  -  PRINT  OUTPUT  OPTION 

-  0,  FOR  FINAL  RESULTS  OUTPUT  ONLY 

-  1,  PRINT  K,M,C  MATRICES  AND  VALUES  OF  ROOTS  AT  VARIOUS 
STAGES  OF  CONVERGENCE 

NOTE:  NODAL  MASS  MATRIX  IS  ADDED  TO  CONSISTENT  MASS  MATRIX  WHICH 
IS  A  FUNCTION  OF  DISTRIBUTED  MASS  DENSITY  (RHO)  INPUT  IN 
EDINPT  (VMP  MATRIX) 


CARD  5  FORMAT  (REQUIRED  IF  ISPIN-1 ) 

SV1.SV2.SV3 
FREE  FORMAT 
WHERE: 

SV1  -  SPIN  RATE  ALONG  GLOBAL  X-AXIS  (RAD/SEC) 
SV2  -  SPIN  RATE  ALONG  GLOBAL  Y-AXIS  (RAD /SEC) 
SV3  -  SPIN  RATE  ALONG  GLOBAL  Z-AXIS  (RAD/SEC) 


CARD  6  FORMAT  (REQUIRED  IF  IPROB  IS  NOT  EQUAL  TO  8) 

I N  DEX ,  N  R ,  INORM ,  PU ,  PL ,  I SOLT ,  IN  DATA 
FREE  FORMAT 
WHERE • 

INDEX  -  INDICATOR  FOR  NUMBER  OF  EIGENVALUES  TO  BE  COMPUTED. 

-  1 ,  WHEN  FIRST  NR  SMALLEST  ROOTS  ARE  NEEDED. 

-  2.  WHEN  ALL  ROOTS  LYING  WITHIN  LIMITS  OF  PU  AND 

PL  ARE  TO  BE  COMPUTED. 

NR  -  NUMBER  OF  EIGENVALUES  TO  BE  COMPUTED. 

INORM  -  THE  PARTICULAR  DEGREE  OF  FREEDOM  AT  WHICH  THE  VECTOR 

SCALOR  IS  USED  TO  NORMALIZE  THE  EIGENVECTORS. 

-  0,  THE  ROUTINE  NORMALIZES  THE  VECTORS  WITH  RESPECT  TO 
THE  ELEMENT  OF  Y  HAVING  LARGEST  MODULUS . 

-  -1,  ELEMENT  Y  OR  YD  IS  USED  FOR  NORMALIZATION. 

PU  -  UPPER  EIGENVALUE  LIMIT 

PL  -  LOWER  EIGENVALUE  LIMIT 

I SOLT  -  VERSION  OF  EIGSOL 

-  0  OR  1 ,  STANDARD  VERSION 

-  2,  SECOND  VERSION 

INDATA  -  INPUT  DATA  OPTION  FOR  EIGSOL 

-  0,  READ  DATA  FROM  FILE 

-  1 ,  READ  DATA  FROM  CARDS 


CARD 


CARD 


CARD 


CARD 


CARD 


7  FORMAT  (REQUIRED  IF  IDRS  «  1) 

TF,NTTS , IUV, DELT, IDDI 
FREE  FORMAT 

TF  -  TOTAL  TIME  SPAN  TO  EVALUATE  DYNAMIC  RESPONSE 

NTTS  -  TOTAL  NUMBER  OF  SETS  OF  LOAD /ACCELERATION  INPUT  DATA 

IUV  -  INDEX  FOR  INITIAL  DISPLACEMENT  AND  VELOCITY  INPUT 

-  0,  NO  INITIAL  DISPLACEMENT/VELOCITY  INPUT 

-  1 ,  IF  EITHER  INITIAL  DISPLACEMENT  OR  VELOCITY  OR  BOTH  ARE 

NON -ZERO 

DELT  -  TIME  INTERVAL  FOR  RESPONSE  CALCULATION 
IDDI  -  INDEX  FOR  DYNAMIC  DATA  INPUT 

-  1,  NODAL  LOAD  INPUT 

-  2,  NODAL  ACCELERATION  INPUT 


8  FORMAT  (REQUIRED  IF  IPROB  *  5  OR  7) 

G 

FREE  FORMAT 
WHERE: 

G  -  STRUCTURAL  DAMPING  PARAMETER  (K=K(1+IG)) 


9  FORMAT  (REQUIRED  IF  IPROB  -  4  OR  5) 
TO  READ  DIAGONAL  VISCOUS  DAMPING  MATRIX 
FORMAT  (6E10.4) 


10  FORMAT  (REQUIRED  IF  IPROB  -  6  OR  7) 
TO  READ  ALPHA, BETA 
FORMAT  (2E10.4) 


11  FORMAT  (REQUIRED  IF  IPROB  -  6  OR  7  AND  ALPHA  -  BETA  -  0.0) 
TO  READ  BANDED  VISCOUS  DAMPING  MATRIX 
FORMAT  (6E10.4) 
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DESCRIPTION  OF  NODAL  DATA  (NODCOR)  INPUT 


(1)  INPUT  NODAL  COORDINATE  DATA  AT  RANDOM 

(2)  THE  FINAL  DATA  IS  AUTOMATICALLY  FORMED  IN  SEQUENCE 

(3)  IN  AUTOMATIC  MESH  GENERATION,  THE  IMPOSED  DISPLACEMENT 
BOUNDARY  CONDITIONS  AND  OTHER  ELEMENT  PROPERTIES  OF 
INTERMEDIATE  NODES  PERTAIN  TO  THE  LAST  CARD  OF  THE 
SEQUENCE. 

(4)  THE  THIRD  POINT  NODES  FOR  LINE  ELEMENTS  LYING  ON 
LOCAL  X-Y  PLANE  OF  AN  ELEMENT  MAY  BE  ANY  ACTIVE  NODE 
OR  DUMMY  ONES  PLACED  AT  THE  END  OF  THE  LIST  WITH 
UX=UY=UZ=TX=TY=TZ=0 . 


CARD  FORMAT 


NODE 

NO. 

NODAL 

X 

COORDINATES 

Y  Z 

NODAL  BOUNDARY 
UX  UY  UZ 

CONDITIONS 

TX  TY  TZ 

INCR 
I  INC 

15 

El  0.4 

El  0.4  El  0.4 

15  15  15 

15  15  15 

15 

NOTE: 

(1)  FOR  NODAL  BOUNDARY  CONDITION 
SET  VALUE  -  0  FOR  FREE, 

-  1  FOR  CONSTRAINED 

(2)  FOR  INCREMENT 

SET  VALUE  -  0  FOR  NO  INCREMENTATION 

-  N  INCREMENT  BY  N  UNTIL  FULL  VALUE 


PARAMETERS  USED  IN  PROGRAM 
(1  )  IN  -  NODE  NUMBER 

(2)  VN (  )  -  NODAL  COORDINATES 

(3)  IIN  (  )  -  NODAL  BOUNDARY  CONDITIONS 

(4)  INC (  )  -  INCREMENT 

(5)  NN  -  TOTAL  NUMBER  OF  NODES 
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EXAMPLE  OF  INPUT: 


10 


INPUT 

1  0.0 

0.0 

0.0 

3  20.0 

0.0 

0.0 

4  0.0 

10.0 

0.0 

10  0.0 

30.0 

0.0 

6  20.0 

10.0 

0.0 

12  20.0 

30.0 

0.0 

5  10.0 

10.0 

0.0 

11  10.0 

30.0 

0.0 

NOTE:  NODES  2,7,9  AND  8  ARE  GENERATED 


oooooooo 


DESCRIPTION  OF  EDINPT  INPUT 


(1)  INPUT  STRUCTURAL  ELEMENT  DATA  AT  RANDOM 
CARD  FORMAT  FOR  SET  (1) 

ELEM  ELEM  NODE  NODE  NODE  NODE  NODE  NODE  NODE  NODE  IMPP  IEPP  ITMPP  IPRR 
TYPE  NO.  1  2345678  /ITHTH 


LINE  *  *  *  **  IEC1  IEC2  FX  + 

1 


SHELL  *****  +  *  ***  *  * 

QUAD 
2 


SHELL  *  *  *  *  +  *  ***  * 

TRIA 
3 


NOTE: 


(1)  **  FOR  THIRD  POINT  (TYPE  1) 

(2)  ***  FOR  ELEMENT  THICKNESS  TYPE  (TYPE  2,3) 

(3)  +  ELEMENT  INCREMENT  NUMBER  FOR  AUTOMATIC  GENERATION  OF 
INTERMEDIATE  ELEMENT  DATA  BY  LINEAR  INTERPOLATION, 

HAVING  PROPERTIES  SAME  AS  LAST  CARD  IN  SEQUENCE 

(4)  FOR  LINE  ELEMENT  -  IEC  -  1  FOR  PIN  ENDED,  -  0  FOR  RIGID  ENDED 

(5)  FX  -  INITIAL  FORCE  IN  LOCAL  X  PLANE  (TYPE  1) 

(6)  IMPP  -  ELEMENT  MATERIAL  PROPERTY  TYPE 

(7)  IEPP  -  ELEMENT  GEOMETRY  PROPERTY  TYPE  (TYPE  1) 

(8)  ITMPP  -  ELEMENT  TEMPERATURE  PROPERTY  TYPE 

(9)  IPRR  -  ELEMENT  PRESSURE  PROPERTY  TYPE 

(10)  ITHTH  -  ELEMENT  THICKNESS  PROPERTY  TYPE  (TYPE  2,3) 


FORMAT  14(15) 


CARD  FORMAT  FOR  SET  (2),  ELEMENT  PROPERTY  CARD  FOR  TYPE  1  ONLY 


PROP 

TYPE 

A  JX 

IY 

IZ 

15 

E10.A  E10.A 

E10.A 

E10.A 

FORMAT 

(I5.3E10.A) 

NOTE: 

(D 

A  -  AREA 

(2)  JX  -  TORSIONAL  MOMENT  OF  INERTIA 

(3)  IY  -  MOMENT  OF  INERTIA 
(A)  IZ  -  MOMENT  OF  INERTIA 


CARD  FORMAT  FOR  SET  (3),  ELEMENT  THICKNESS  DATA  CARD  FOR  TYPE  2,3  ONLY 

THICK  T 
TYPE 

15  E10.A 

FORMAT  (I5.E10.A) 

NOTE: 

(1)  T  -  THICKNESS 


CARD  FORMAT  FOR  SET  (A),  ELEMENT  MATERIAL  PROPERTY  CARD 
[TO  BE  SET  FOR  GENERAL  MATERIAL  CASE] 


MATL  ELASTIC  CONSTANTS 
TYPE 

FOR  ISOTROPIC 


* 

*<E) 

*(MU) 

* (ALP ) 

*(RHO) 

FOR 

ORTHOTROPIC 

• 

* 

* 

* 

* 

* 

* 

* 

* 

* 

* 

*(ALPX) 

*(ALPY) 

*(RHO) 

FOR 

ANISOTROPIC 

* 

* 

* 

* 

★ 

* 

★ 

* 

* 

* 

* 

* 

* 

* 

★ 

* 

* 

* 

★ 

* 

* 

★ 

*(ALPX) 

* (ALPY ) 

*(ALPXY) 

*(RHO) 

FORMAT  (15, 7E1 0. A/5X, 7E1 0. A/5X, 7E1 0. A/5X, AE1 0. A) 
NOTE:  DATA  INPUT  IS  REQUIRED  FOR  EACH  MATERIAL  CASE 


CARD  FORMAT  FOR  SET  (5),  ELEMENT  TEMPERATURE  INCREMENT  CARD 

TEMP  T  DTDY  DTDZ  TEMP  T  DTDY  DTDZ 

TYPE  TYPE 

15  El  0.4  El  0.4  El  0.4  15  E10.4  E10.4  E10.4 

FORMAT  2(I5,3E10.4) 

NOTE: 

(1)  T  -  TEMPERATURE  INCREMENT 

(2)  DTDY  -  Y-TEMPERATURE  GRADIENT 

(3)  DTDZ  -  Z-TEMPERATURE  GRADIENT 


CARD  FORMAT  FOR  SET  (6),  ELEMENT  PRESSURE  LOAD  CARD 

PRESS  PRESS  PRESS  PRESS  PRESS  PRESS  PRESS  PRESS  PRESS  PRESS 

TYPE  TYPE  TYPE  TYPE  TYPE 

15  El  0.4  15  El  0.4  15  E10.4  15  E10.4  15  E10.4 

FORMAT  5(15, E10. 4) 

NOTE: 

ELEMENT  TYPE  1  -  UNIFORM  PRESSURE  ALLOWED  IN 

LOCAL  Y-DIRECTION  ONLY 

ELEMENT  TYPE  2,3  -  UNIFORM  PRESSURE  ALLOWED  IN 

LOCAL  Z-DIRECTION  ONLY 


DESCRIPTION  OF  GCINPT  INPUT 


SET  (1  )  INPUT  NODAL  FORCE  (P)  OR  MASS  (M)  MATRIX  AT  RANDOM 
(ONLY  IF  IPLUMP /IMLUMP  NOT  EQUAL  TO  0) 

SET  (2)  INPUT  DYNAMIC  LOADING  (ONLY  IF  IDRS  -  1) 

SET  (3)  INPUT  NODAL  CONNECTIVITY  AT  RANDOM 


CARD  FORMAT  FOR  SET  (1) 

NODE  NUMBER  DOF  P  OR  M 

15  15  El  0.4 

FORMAT  (215, El 0.4) 

NOTE:  (1)  FOR  NODAL  LOAD,  EACH  CASE  IS  TERMINATED  BY  SETTING  THE 

NODE  NUMBER  FOR  THE  NEXT  CASE  TO  A  NEGATIVE  NUMBER,  SAY  -1 
(2)  FOR  IPROB  -  8,  IF  IPLUMP  -  0  THEN  TERMINATE  INPUT  WITH  -1 


CARD  FORMAT  FOR  SET  (2) 

CARD  1  FORMAT  (REQUIRED  IF  IUV  -  1  AND  IDRS  -  1) 

TO  READ  INITIAL  DISPLACEMENT/VELOCITY  AT  RANDOM,  TERMINATED  BY  -1 
FORMAT  (2I5.2E15.5) 

NODE  NUMBER  DOF  U(I)  V(l) 

*  *  *  * 


.  ... 

*  *  *  * 

-1 

CARD  2  FORMAT  (REQUIRED  IF  IDRS  -  1  AND  NTTS  >  0) 

TO  READ  NTTS  NUMBER  OF  SETS  OF  NODAL  LOAD /ACCELERATION  DATA) 
FORMAT  (El  5. 5) 

TZ 

* 

FORMAT  (215, El 5.5) 

NODE  NUMBER  DOF  LOAD /ACCELERATION 

★  *  h 


-1 

NOTE:  REPEAT  CARD  2  DATA  AS  AVOVE  FOR  NTTS  NUMBER  OF  SETS 
EACH  TERMINATED  BY  -1 


CARD  FORMAT  FOR  SET  (3) 

4(14,11 ,14,11  ,E10. 4) 

(NODE  DOF  NODE  DOF  CONNECTIVITY)  4  SETS  PER  ROW 

COEFFICIENT 


3-12 


FURTHER  MAIN I  INPUT 


INPUT  OF  VISCOUS  DAMPING  MATRICES 


CARD  1  FORMAT  (REQUIRED  IF  IPROB-4  OR  5) 

TO  READ  DIAGONAL  VISCOUS  DAMPING  MATRIX , C (N , 1 )  IN  GCS 
FORMAT:  (6E10.4) 


CARD  2  FORMAT  (REQUIRED  IF  IPROB-6  OR  7) 
TO  READ  'ALPHA'  AND  'BETA' 

SO  THAT  [C]  -  ALPHA* [K]  +  BETA*[M] 
FORMAT:  (2E10.4) 


CARD  3  FORMAT  (REQUIRED  IF  IPROB-6  OR  7  AND  ALPHA-BETA-0.0) 
TO  READ  BANDED  VISCOUS  DAMPING  MATRIX, C(N , Ml  1 )  IN  GCS 
FORMAT:  (6E10.4) 


NOTE:  DATA  IN  CARDS  1  AND  3  MUST  CONFORM  TO  'ZDBC.FDBC  AND  IDBC' 
INHERENT  IN  PROBLEM  UNDER  CONSIDERATION. 

WHERE  * 

ZDBC  -  ZERO  DEFLECTION  BOUNDARY  CONDITIONS,  INPUT  IN  NODCOR 
FDBE  -  FINITE  DEFLECTION  BOUNDARY  CONDITIONS, 

INPUT' IN  SET  (3)  IN  GCINPT 

IDBE  -  INTERDEPENDENT  DEFLECTION  BOUNDARY  CONDITIONS, 

INPUT  IN  SET  (3)  IN  GCINPT 


PARAMETERS  USED  IN  THE  PROGRAM 


(1)  NEL 

(2)  NEP 

(3)  NET 

(4)  NMAT 

(5)  NMECN 


(6)  NTMP 

(7)  NPR 

(8)  IEP (  )  - 

(9)  VEP (  ,  )  - 

(10)  IET(  ) 

(11)  VET(  ) 

(12)  IMP (  )  « 

(13)  VMP (  ,  )  - 

(14)  ITMP(  )  - 

(15)  VTMP (  ,  )- 

(16)  I PR (  )  - 

(17)  VPR(  )  - 


TOTAL  NUMBER  OF  ELEMENTS 

TOTAL  NUMBER  OF  ELEMENT  PROPERTY  TYPES  (LINE  ELEMENTS) 

TOTAL  NUMBER  OF  ELEMENT  THICKNESS  TYPES  (SHELL  ELEMENTS) 

TOTAL  NUMBER  OF  ELEMENT  MATERIAL  TYPES 

NUMBER  OF  MATERIAL  ELASTIC  CONSTANTS 

4,12,25  FOR  ELASTIC  ISOTROPIC,  ORTHOTROPIC,  OR 

ANISOTROPIC  CASE,  RESPECTIVELY 

TOTAL  NUMBER  OF  ELEMENT  TEMPERATURE  TYPES 

TOTAL  NUMBER  OF  ELEMENT  PRESSURE  TYPES 

ELEMENT  PROPERTY  TYPE  NUMBER,  LINE  ELEMENTS  (TYPE  1) 

ELEMENT  PROPERTIES 

ELEMENT  THICKNESS  TYPE  NUMBER,  SHELL  ELEMENTS  (TYPE  2,3) 

ELEMENT  THICKNESS 

ELEMENT  MATERIAL  TYPE  NUMBER 

ELEMENT  ELASTIC  CONSTANTS 

ELEMENT  TEMPERATURE  TYPE  NUMBER 

ELEMENT  TEMPERATURE  DATA 

ELEMENT  PRESSURE  TYPE  NUMBER 

ELEMENT  PRESSURE  DATA 
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DEVELOPMENT  OF  A  UNIFIED  NUMERICAL 
PROCEDURE  FOR  FREE  VIBRATION  ANALYSIS 

OF  STRUCTURES'!* 

K.  K.  GUPTA* 

Jet  Propulsion  Laboratory,  California  Institute  of  Technology,  Pasadena,  California,  U.S.A. 

SUMMARY 

This  paper  presents  the  details  of  a  unified  numerical  algorithm  and  the  associated  computer  program 
developed  for  the  efficient  determination  of  natural  frequencies  and  modes  of  free  vibration  of  structures. 
Both  spinning  and  nonspinning  structures  with  or  without  viscous  and/or  structural  damping  may  be  solved 
by  the  current  procedure;  in  addition,  the  program  is  capable  of  solving  static  problems  with  multiple  load 
cases  as  well  as  the  quadratic  matrix  eigenproblem  associated  with  a  finite  dynamic  element  structural 
discretization.  A  special  symmetric  matrix  decomposition  scheme  has  been  adopted  for  matrix  tri- 
angularization,  which  renders  the  program  rather  efficient  and  economical.  Also,  a  novel  bisection  scheme 
is  described  that  further  accelerates  the  solution  convergence  rate,  particularly  for  the  case  of  repeated 
roots. 

The  associated  computer  program  adopts  an  out-of-core  solution  strategy,  thereby  enabling  effective 
solutions  to  be  achieved  for  large,  complex,  practical  structures.  A  complete  listing  of  the  program  written 
in  FORTRAN  V,  for  the  UNIVAC  1100/82  computer,  along  with  the  source  deck  is  available  for  ready 
use. 

INTRODUCTION 

The  dynamic  response  analysis  is  of  primary  importance  in  the  design  of  a  wide  range  of  practical 
structures,  such  as  spacecraft,  buildings,  and  rotating  machineries,  among  others.  A  vital 
preliminary  for  such  an  analysis  involves  the  determination  of  the  natural  frequencies  and  the 
associated  modes.  This  is  achieved,  first,  by  discretizing  the  continuum  by  a  standard  technique, 
such  as  the  finite  element  method,  yielding  simultaneous  algebraic  equations;  the  resulting 
eigenvalue  problem  is  then  suitably  solved  to  yield  the  desired  roots  and  vectors.  For  most 
complex  practical  structures,  such  an  idealization  results  in  a  rather  large  number  of  simul¬ 
taneous  equations,  which  are  usually  of  highly  banded  configurations.  In  order  to  effect  an 
economical  solution,  the  associated  eigenproblem  analysis  routine  must  be  designed  to  fully 
exploit  the  inherent  matrix  sparsity.  Furthermore,  due  to  the  limited  core  storage  available  in 
present  computers,  it  is  advantageous  to  adopt  an  out-of-core  solution  strategy  that  provides 
effective  solutions  for  practical  structures  of  almost  any  magnitude  add  complexity. 

While  many  structures  are  nonrotating  in  nature,  some  are  subjected  to  uniform  rotations. 
Abo,  such  structures  may  exhibit  viscous  or  structural  damping  or  a  combination  of  both.  The 
associated  eigenvalue  problems  are  characterized  by  distinctive  matrix  equations.  Furthermore, 
when  finite  dynamic  elements  are  used  for  structural  discretization,  a  quadratic  matrix 

tThe  research  described  in  this  paper  was  carried  out  by  the  Jet  Propulsion  Laboratory,  California  Institute  of 
Technology,  and  was  sponsored  jointly  by  the  Air  Force  Office  of  Scientific  Research  (AFOSR)  and  the  Large  Space 
Structures  Technology  (LSST)  Project  Office  at  the  NASA  Langiey  Research  Center. 
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eigenvalue  formulation  is  involved.  A  general  formulation  unifying  the  various  eigenvalue 
problems  may  be  presented  as 

M$+C4+Kq-0  (1) 

in  which  M,  C  and  K  are,  in  general,  the  inertia,  damping  and  stiffness  matrices,  respectively. 
The  individual  eigenproblems  are  identified  next. 

Case  /.  Undamped  free  vibration  and  buckling 

The  related  matrix  formulation  becomes 

Mq+Kq-0  (2) 

the  solution  of  which  is  taken  to  be  q  -  eu';  equation  (2)  then  reduces  to 

(Ke-A2M)q-0  (2a) 

where  KE  is  the  elastic  stiffness  matrix,  and  A  the  natural  frequencies.  The  associated  buckling 
problem  is  characterized  by 

(KE-AKo)q-Q  (2b) 

in  which  Ka  is  the  geometrical  stiffness  matrix  and  A  the  compressive  buckling  load.  The 
formulation  for  the  associated  problem  of  tree  vibration  of  prestressed  structures  is  given  by 

[(KB-Ko)-A2M]q-0  (2c) 

in  which  the  compressive  load  is  assumed  to  have  a  positive  sign.  Equations  (2aH2c)  are 
characterized  by  real  roots  and  vectors  since  Kb,  K <*  M  are  real,  symmetric  matrices,  KE  also 
being  positive  definite  in  nature;  for  structures  exhibiting  rigid  body  motion,  a  non-negative 
definite  form  of  Ke  is  obtained. 

Case  II.  Undamped  free  vibration  of  spinning  structures 

For  undamped  structures  spinning  at  a  uniform  rate  fl,  equation  (1)  assumes  the  form 

M4+Ce4+Kq-0  (3) 

where  C«  »  skew -symmetric  Coriolis  matrix,  being  a  function  of  ft;  K  »  KE +K0  +  S';  Ko  and 
K'  are  the  geometrical  stiffness  and  centrifugal  forces  matrices,  respectively,  both  being 
functions  of  fl2.  The  solution  (3)  is  assumed  as  q  -  e*,  and  the  resulting  eigenvalue  problem 
takes  the  form1 

(K+pCe+p2M)q-0  (3a) 

in  which  p  is  pure  imaginary,  such  roots  and  associated  complex  vectors  occurring  in  conjugate 
pairs.  Further,  K  and  M  are  assumed  to  be  symmetric  and  positive  definite  for  small  vibrations. 

Case  m.  Quadratic  matrix  equations  ■ 

If  structural  discretization  is  achieved  by  finite  dynamic  elements  (FDEs),  the  resulting 
frequency-dependent  stiffness  and  inertia  matrices  are,  first,  expressed  in  terms  of  ascending 
powers  of  the  frequencies  A : 


K-Ko  +  A4K«  + 
M-Mo  +  A2M,+ 


(4) 
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resulting  in  the  quadratic  matrix  eigenvalue  problem 

[Ko  -  A  2Mo  -  A  4(M2  -  K«)]q  -  0 


(4a) 


where  Ko»K e,  Mo  *  M  are  the  usual  elastic  stiffness  and  mass  matrices;  M2,  Kt  are  the 
higher  order  dynamic  corrections  terms,  the  inclusion  of  which  in  the  free  vibration  analysis  is 
known  to  effect  a  significant  increase  in  the  solution  convergence  rate.2  Equation  (4a)  is  of 
similar  form  to  equation  (3a),  but  possesses  real  roots  and  vectors.  Equation  (4a)  may  also  be 
written  as 


(K-A2M-A4£)q-0 

which  may  be  further  rearranged  as 


with  q=«  A2q  and  which  is  of  the  form 

(E-A2F)y«0 

where  E  and  F  are  symmetric  matrices,  E  also  being  positive  definite  in  nature. 


(4b) 

(4c) 

(4d) 


Cases  IV  and  V.  Damped  free  vibration  of  spinning  structures 

The  associated  equations  of  free  vibration  without  or  with  structural  damping  are  expressed  as 

Mlj+tC.+CrfN+Kaq-O  (5) 

Mq+(Cc  +  C-)4+KH(l  +  «*g)q-0  (6) 

for  cases  IV  and  V,  respectively,  in  which  Cd  is  the  viscous  damping  matrix  assumed  to  be 
diagonal,  i*  is  the  imaginary  number  V-l,  and  g  is  the  structural  damping  parameter. 
Substituting  q  -  e*  in  the  above  equations,  the  resulting  eigenvalue  problems  have  the  following 
form: 

[Ke+p(Ce  +  C,i)+p2M]q-0  (5a) 

[K«(l+i*g)+p(Ce+C,)VM]q-0  (6a) 

in  which  the  roots  p,  as  well  as  the  associated  vectors,  are  obtained  as  complex  conjugate  pairs. 

Cases  VI  and  VII.  Damped  free  vibration  of  nonspinning  structures 
Equation  (1)  assumes  the  following  form  for  cases  VI  and  VII,  respectively: 


M<i  +  C4  +  Keq-0  (7) 

M<i+ Cq + Kg(l  +  /*g)q  ■  0  (8) 

and  the  related  eigenproblem  formulations  are  defined  as 

[KH+pC+p2M]q-0  (7a) 

[Kg(  1  +  i'*g) + pC + p2M]q  ■  0  (8a) 


where  the  roots  and  vectors  have  forms  similar  to  those  pertaining  to  cases  IV  and  V;  C  is  the 
viscous  damping  matrix  of  general  banded  form. 
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Case  Vm.  Simultaneous  equations 
Solution  of  simultaneous  equations 

Sq-b  (9) 

is  effected,  where  S  is  either  real,  symmetric  or  Hermitian  in  nature,  and  b  is  a  set  of  arbitrary 
vectors.  An  effective  solution  of  equation  (9)  is  essential  in  the  analysis  of  the  various  eigenvalue 
problems. 

The  purpose  of  this  paper  is  to  present  the  details  of  a  unified  numerical  algorithm  and  an 
associated  computer  program  developed  for  the  solution  of  the  above  eigenvalue  problems 
defined  by  cases  I  to  VII.  Such  a  program  fully  exploits  the  banded  nature  of  the  related  matrices 
and  enables  computation  of  only  a  few  desired  roots  lying  within  a  specified  bound  without 
having  to  compute  any  other.  Since  the  computer  program  employs  an  out-of-core  solution 
strategy,  it  enables  effective  solutions  to  be  obtained  for  complex  practical  problems  of  rather 
large  magnitudes.  Numerical  results  are  also  presented  in  detail,  testifying  to  the  relative 
efficiency  of  the  present  procedure.  This  is  followed  by  a  summary  of  conclusions. 

BASIC  NUMERICAL  SCHEMES 

A  solution  of  the  general  eigenvalue  problem  defined  by  equation  (1)  is  obtained  by  first 
rearranging  the  matrices  as  follows: 


which  may  also  be  written  as 

By+Ay-0  (11) 

where 


Substituting  y  «  e*  as  its  solution,  equation  (11)  takes  the  form 

(B-t-pA)y-O  (12) 

which  may  also  be  written  as 

(B-A  A*)y-0  (13) 

where  A  « i*p  and  A*  -  i* A. 

Isolation  of  roots 

For  the  particular  case  when  A*  is  a  Hermitian  matrix  and  B  is  real,  symmetric  and  positive  or 
non-negative  definite  in  nature,  the  Sturm  sequence  property1  is  valid  for  the  formulation 
depicted  by  equation  ( 13),  the  associated  roots  being  real  in  nature.  Thus,  for  a  given  value  of  A, 
the  number  of  changes  in  the  signs  of  the  leading  principal  minors  f,( A )  is  equal  to  the  number  of 
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roots  ofB-AA*  having  algebraic  values  less  than  A.  As  a  special  case,  this  property  is  also  valid 
for  the  formulation  involving  a  pair  of  real  symmetric  matrices.  Such  a  property  enables  the 
isolation  of  only  a  few  desired  roots  lying  within  a  specified  upper  and  lower  bound  [AM,  A »],  by  a 
repeated  bisection  technique,  without  having  to  compute  any  other.  Thus,  at  a  particular  stage  of 


A 


J 


Figure  1.  Convergence  scheme  for  repeated  roots 


computation,  the  Sturm  count  determines  the  number  of  roots  of  B- A  A*,  say,  x,  lying  within 
[Am  Af],  with  Am  « (Au  +  a()/2.  Ah,  A<  being  the  current  upper  and  lower  bounds,  respectively. 
The  upper  bound  of  the  x  roots  and  the  lower  bounds  of  the  rest  are  then  set  to  Am,  the  latter 
being  implemented  only  if  their  current  lower  bounds  are  smaller  than  Am.  The  desired  roots  are 
automatically  isolated  and  their  individual  bounds  determined  when  this  process  is  continued. 
Furthermore,  if  a  number  of  roots  are  found  within  the  bound  [Am,  A(],  such  that  the  absolute 
value  of  (Am  -  Ai)/A|  is  less  than  the  root  separation  parameter  EPS,  then  they  are  considered  as 
repeated  ones  with  a  numerical  value  equal  to  Am.  However,  the  latter  process  tends  to  be  rather 
slow  for  extracting  repeated  or  close  roots.  A  novel  strategy  has  been  developed  in  connection 
with  the  present  work  which  essentially  reduces  the  root  extraction  time  for  repeated  roots  to 
that  of  distinct  ones.  Thus,  during  the  bisection  process,  when  monitoring  the  bounds  of  a  group 
of  roots,  if  the  number  of  such  roots,  say,  r,  remain  the  same  while  their  upper  and  lower  bounds 
each  changes  at  least  once,  then  a  multiplicity  test  is  immediately  carried  out.  The  inverse 
iteration  process  described  in  the  next  section  is  then  used,  employing  the  respective  upper  and 
lower  bound  values  to  accurately  locate  two  roots  by  converging  from  both  ends.  If  these  two 
root  values  are  found  to  be  identical,  the  r  roots  are  then  asumed  to  be  multiple  ones  having  the 
numerical  value  as  that  of  the  two  converged  roots.  If,  on  the  other  hand,  the  two  converged 
roots  are  distinct  in  values,  they  are  accepted  as  the  true  values  of  the  respective  roots  and  the 
bisection  process  is  continued  as  usual.  This  special  procedure  is  depicted  in  Figure  1,  where  the 
roots  Ak  to  A„  are  repeated  in  nature.  First  they  are  isolated  within  bounds  [A*  Aj]  when  two 
more  bisections  are  needed  to  satisfy  the  criterion  for  the  present  strategy.  The  final  bounds  A  'k 
and  A  *  are  then  utilized  to  converge  from  both  sides.  This  current  procedure  has  been  found  to 
be  much  superior  than  the  usual  repeated  bisection  technique  for  multiple  roots. 

For  the  present  set  of  problems,  the  Sturm  sequence  count  involving  the  number  of  changes  in 
the  signs  of  the  leading  principal  minors  is  equivalent  to  counting  the  number  of  negative 
diagonal  elements  of  the  decomposed  matrix.  Depending  on  the  type  of  problem,  such  an 
operation  is  performed  on  the  following  matrix  formulation: 


Cases  /,  VI  and  VII 

Matrix  triangularization  is  performed  on 

K-A:M 


when  all  operations  involve  real  numbers. 


(14) 
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V 
* ; 
-  \ 


v'; 


h 

V 


Casein 

Similar  matrix  manipulations  are  executed  with  the  matrix 

K-A2M-A4£  (15) 

also  involving  real  numbers  only. 

Cases  II,  IV  and  V 
The  decomposition  of  the  matrix 

K-/*ACc  — A2M  (16) 

is  implemented,  when  all  operations  are  performed  in  complex  arithmetic. 

Location  of  roots  and  computation  of  vectors 

Once  the  roots  are  isolated  by  the  repeated  bisection  procedure,  an  inverse  iteration 
technique  is  adopted  for  the  simultaneous  determination  of  individual  roots  and  vectors.3'4  In 
this  process  the  middle  point  of  the  bounds  of  the  isolated  rth  root  is  taken  as  the  starting  root 
iteration  value: 

A'„»(AL+Aj)/2  (17) 

which  is  utilized  to  effect  the  triangularization  of  the  relevant  left-hand  side  matrix.  A  starting 
vector  is  then  chosen  to  consist  entirely  of  unit  real  scalars  to  start  matrix  iterations;  for  complex 
operations,  the  imaginary  parts  of  the  scalars  are  assumed  to  be  zero.  At  the  end  of  each  iteration 
a  Rayleigh  quotient  is  used  to  obtain  a  new  estimate  of  the  root  under  consideration.  The  matrix 
formulations  adopted  for  the  inverse  iteration  procedure  for  a  typical  iteration  are  summarized 
next: 

Case  I 

[K  -  (A  (18) 

Ni+i  being  a  normalizing  factor. 

Case  II 

Lower  half  of  equation  (13): 

[K - i*\rmCe  -  (A ^)2M]q;+l  -  (V(+l[A ;Mq[+  /*Mq,'+  i*Ccq'(]  (19) 

Upper  half: 

[M^l+i  +  i*A mMql+i]  “  ~N<+ 1 1'Mq,'  (19a) 


Case  m 

Lower  half  of  equation  (4d): 

[K-(A,m)2M-(A;)4e]q;+I  MqJ+iA^fcqn  (20) 


Upper  half: 


6^i-(Ai.)2dq('*,-Af,+,e:q; 


(20a) 
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Cases  TV  and  V 
Lower  half  of  equation  (13): 

[K — i*A  i,Ce — (A  !*)2!V«]qJ+i  “  N)+i[A  *Mqi' + i*Mqj + i*Cqi]  (21) 

in  which  C  »  Cc  +  Cd. 

Upper  half:  As  in  case  II 

Cases  M  and  MI 
Lower  half  of  equation  (13): 

[K — (A  *)2M]q[+i »  Ni+i[A  mMq<  +  (*M4i + /*Cqj]  (22) 

in  which  C  is  the  viscous  damping  matrix  of  usual  banded  form. 

Upper  half:  As  in  case  II. 

For  each  root,  triangularization  of  the  left-hand  side  of  the  above  equations  is  performed  only 
once  at  the  beginning  of  the  iteration.  In  subsequent  iteration  steps,  their  solutions  are  achieved 
by  the  simple  back-substitution  process. 

The  Rayleigh  quotient  is  used  at  each  iteration  step  to  achieve  a  new  estimate  for  the  root; 
their  detailed  expressions  are  presented  below: 


Case  I 


(Ai+i  )2  *  (qI+i)TKqi+i/(q<+i)TMqi+i 

(23) 

Case  II 

»(tf+t)TByJ*,/(tf+i)TA*yU. 

(24) 

where 

(ymfByl+t  -  (i!+i)TM*U  +  (tf*i)TKqU 

(24a) 

(yJ+l)TA*y[+i  -  -  »*( +  i*(tf +,)TMtf*i  +  «*(fcV  i)TC.q,r+i 

(24b) 

Casein 

(a;+1  )2  -  (y;+,)TEy;*,/(y;+,)TFy;+l 

(25) 

where 

(y<+i)TEyj+i  ■  (f|i+i)T^i+i  +  (q«+i)TKq«+i 

(25a) 

(y«.i)TlW*i  “  (q<+i)TCf|(+i  +(f|(+i)T6qJ+i  +(q<+i)TMqj+i 

(25b) 

Case  W 

As  in  case  II,  but  C«  is  replaced  by  €e  +  Cd. 
Case  V 

As  in  case  IV,  but  K  is  replaced  by  K(1  +  i*g). 
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Case  VI 

As  in  case  II  with  C  replacing  C* 

Case  VII 

As  in  case  VI,  but  K  is  changed  to  K(1  +  i*g). 

The  new  estimate  for  A /+i,  obtained  by  using  the  Rayleigh  quotient,  is  next  used  to  check  the 
pattern  of  root  convergence.  Thus,  if  the  factor 

A  |A|  +  1  A, |  /AAS\ 

ACC »  .  .  (26) 

|A(+ll 

is  found  to  be  smaller  than  the  specified  root  accuracy  parameters  EPS1,  A  <+i  is  then  accepted  as 
the  true  value  of  A';  otherwise,  the  above  steps  involving  inverse  iterations  and  calculation  of 
Rayleigh  quotients  are  continued  until  adequate  convergence  is  achieved. 

Repeated  real  roots  are  determined  entirely  by  the  novel  bisection  process  described  earlier. 
To  achieve  the  associated  vectors,  these  roots  are  first  artificially  separated  relative  to  their 
respective  values  by  an  amount  3  x  2~(/3,  t  being  the  number  of  digits  after  the  binary  point  The 
inverse  iteration  procedure  is  then  directly  applied  to  yield  the  vectors,  when,  at  most,  two 
iterations  are  required  for  their  accurate  determination.  Due  to  extreme  sensitivity  of  eigen¬ 
vectors  to  small  perturbations  in  the  vicinity  of  multiple  roots,  this  procedure  yields  sets  of 
independent  vectors  corresponding  to  such  roots.  The  standard  Schmidt  orthogonalization 
technique  may  then  be  applied  to  transform  the  independent  vectors  into  orthogonal  ones. 


Numerical  stability 

The  triangularization  of  equations  ( 14H22)  has  been  achieved  by  a  symmetric  decomposition 
scheme  that  omits  row  interchanges,  which,  in  turn,  preserves  the  bandwidth  of  the  relevant 
matrices.  This  reduces  the  imputation  time  for  roots  and  vectors  by  a  factor  between  2  and  3, 
when  compared  with  similar  routines  developed  earlier.3'4  Since  the  matrix  equations  (14H22) 
are,  in  general,  not  positive  definite  in  nature,  suitable  pivoting  was  considered  to  be  necessary 
to  preserve  numerical  stability.  In  the  present  work,  some  alternative  measures  have  been 
implemented  in  the  program  that  proves  to  be  effective  in  preserving  numerical  stability. 

Thus,  provided  the  roots  are  not  required  to  be  calculated  to  a  high  precision,  the  resulting 
bisection  process  is  remarkably  stable  in  nature  and  the  chances  of  a  breakdown  are  rather 
slight.3  The  root  accuracy  parameters  EPS  and  EPS1  are  thus  set  to  0-0001  and  0-001, 
respectively,  which  will  result  in  eigenvalues  that  are  sufficiently  accurate  for  problems  usually 
encountered  in  practice.  Also,  the  program  may  easily  be  run  in  double  precision,  if  desired,  so 
that  a  large  number  of  significant  figures  are  retained  in  subsequent  computations.  Moreover,  if 
a  zero  pivot  is  encountered  on  any  rare  occasion,  the  matrix  triangularization  may  be  repeated 
with  a  slightly  perturbed  value  of  the  current  A;  however,  since  the  bisection  process  has 
proved  to  be  rather  well-conditioned,  these  extra  precautions  are  deemed  to  be  unnecessary. 

During  the  inverse  iteration  procedure,  if  a  zero  pivot  is  encountered,  the  program  automa¬ 
tically  replaces  it  by  e,  the  normal  rounding  error,  and  continues  the  process  onward.  A  number 
of  test  cases  have  been  run  using  the  program  and  their  convergence  characteristics  carefully 
monitored.  These  results  indicate  that  the  program  is  accurate  and  reliable  in  nature. 
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DEVELOPMENT  OF  A  COMPUTER  PROGRAM 

A  computer  program  EIGSOLt  (EIGen  problem  SOLution  routine),  based  on  the  numerical 
formulations  presented  in  the  earlier  section,  has  been  developed  to  achieve  efficient  solutions 
to  structural  free  vibration  problems.  Both  spinning  and  nonspinning  structures,  with  or  without 
viscous  and  structural  damping,  as  well  as  the  quadratic  matrix  eigenproblem  associated  with  a 
finite  dynamic  element  discretization,  may  be  analysed  by  the  program.  For  spinning  structures 
the  program  is  limited  to  diagonal  viscous  damping  matrices.  An  adequate  number  of  comment 
cards,  included  in  the  listing,  renders  the  program  to  be  self-explanatory  in  nature. 

An  important  aspect  of  the  analysis  is  related  to  the  solution  of  simultaneous  equations,  either 
real,  symmetric  or  Hermidan, ,  using  an  out-of-core  solution  strategy,  and  this  is  effected  by  the 
subroutine  BANMAT  (BANded  MATrix  solution).  The  relevant  data  are  all  stored  in  secon¬ 
dary  storage  units,  such  as  discs,  which  are  brought  in  the  core  in  suitable  predetermined  block 
format.  A  minimum  core  storage  requirement  of  (Ml  1,  Ml  1)  and  Ml  1  pertaining  to  D  and  AD 
matrices  is  required  to  operate  the  program  BANMAT;  Mil  is  the  half  bandwidth  of  K, 
including  the  diagonal.  The  program  is  designed  to  run  on  the  260KUNTVAC  1100/82 
computer,  which  allows  usage  of  up  to  about  175K  core  storage,  that  enables  achieving  solutions 
to  practical  problems  of  rather  large  magnitude. 

The  main  driver  routine  EIGSOL  repeatedly  calls  subroutine  INPUT  to  effect  data  input 
Thus,  the  upper  symmetric  halves  of  K(N,  Ml  1),  M(N,  MB),  C  (or  Cc)  (N,  MC)  and  C(N,  MCD) 
are  read  in  predetermined  block  format;  N  is  the  order  of  the  matrices,  MB,  MC,  MCD  being  the 
half-bandwidth  of  M,  C,  Cd,  respectively.  Figure  2  shows  the  arrangement  of  data  blocks,  the 


(a)  INPUT  DATA  STORAGE  (b)  SOLUTION  DATA  OUTPUT  STORAGE 

Figure  2.  Data  block  set-up  (or  K.  M  and  C  matrices 


number  of  such  data  blocks  (NDBLK)  being  dependent  on  the  available  core  area  specified  by 
the  parameter  NAC,  which  is  to  be  provided  by  the  user.  Besides  the  main  store  D  of  dimensions 
(Ml  1,  Ml  1),  each  block  will  have  NRAD  number  of  rows,  the  number  of  such  rows  in  the  last 
data  block  being  defined  by  NRLDBK.  Once  the  solution  has  been  achieved,  it  is  then  stored  in 
the  block  format,  as  shown  in  Figure  2(b). 

t  The  physical  program  EIGSOL  is  available  from  the  Computer  Management  and  Information  Center  (COSMIC),  the 
NASA  agency  for  distribution  of  computer  programs. 


Isolation  of  the  desired  roots  is  achieved  by  the  repeated  bisection  technique  effected  by  the 
subroutine  BISECN,  which,  in  turn,  repeatedly  calls  the  subroutine  EIGNV  to  determine  the 
number  of  roots  of  the  present  problem  having  algebraic  values  less  than  the  current  root  value 
under  consideration.  The  subroutine  BANMAT  is  used  by  EIGNV  for  the  Sturm  sequence 
count  of  roots.  Once  the  roots  are  isolated,  they  are  accurately  located  by  the  subroutine 
VECTOR  employing  the  inverse  iteration  procedure,  which  also  simultaneously  yields  the 
associated  vectors.  Two  other  subroutines,  MULT  and  VMULT,  are  also  used  by  the  program 
for  appropriate  matrix  and  vector  multiplications. 
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Figure  3.  Amngment  of  date  in  common  block  my  A  (variables  defined  in  program  listing) 


Figure  3  depicts  the  schematic  arrangement  of  data  in  the  main  common  block  core, 
containing  an  array  A  that  contains  all  major  vectors  and  matrices.  The  starting  addresses  of 
theses  arrays  relative  to  the  array  A  are  also  shown  in  the  figure,  which  are  used  in  the  arguments 
of  the  various  subroutines  called  by  the  main  driver  routine  to  effect  appropriate  equivalence  of 
these  arrays  with  array  A.  This  common  block  is  designed  in  such  a  way  that  it  occupies  the  last 
portion  of  the  computer  data  bank  (DB  ANK).  Thus,  as  long  as  A  has  at  least  a  starting  address  in 
the  first  core  module  of  the  computer,  it  will  automatically  spill  over  to  the  other  data  banks, 
enabling  utilization  of  about  175K  core  memory  in  a  260K  UNIVAC 1 100/82  machine.  Thus, 
the  program  is  capable  of  solving  rather  large  practical  problems.  Depending  on  the  nature  of 
the  problem,  the  program  operates  either  in  real  or  complex  mode. 


NUMERICAL  RESULTS 

An  extensive  number  of  test  cases  have  been  solved  by  the  program  EIGSOL  to  check  out  the 
various  capabilities  offered  by  the  procedure  as  well  as  to  establish  the  relative  efficacy  of  the 
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program  in  comparison  to  other  existing  ones.  Thus,  the  spinning  cantilever  beam  problem, 
presented  earlier,3  is  again  chosen  as  a  numerical  example.  The  basic  elastic  properties  of  the 
beam,  divided  in  10  discrete  elements  of  length  /  and  expressed  in  the  inch-pound-second  unit 
system,  are  as  follows: 

Moment  erf  inertia  (Y-axis) 

Moment  of  inertia  (Z-axis)  -  & 

Area  of  cross-section  - 1-0 

Young’s  modulus  ■  30  x  106 

Nodal  mass  in  translation  *  1 

Nodal  mass  moment  of  inertia 
Scalar  viscous  damping  ■  0*628318 

Structure  damping  parameter  (g)  -  0*01 

Element  length  (/)  *  6 

where  the  direction  of  the  Y-axis  is  chosen  along  the  length  of  the  beam.  A  number  of  computer 
runs  were  performed  corresponding  to  problem  cases  H,  IV  and  V  by  setting  the  appropriate 
input  parameter  IPROB  to  2, 4  and  S,  respectively.  The  beam  was  subjected  to  a  uniform  spin 
rate  0*0-1  Hz,  along  the  Y-axis,  at  the  built-in  end;  the  results  of  such  analyses  are 
summarized  in  Table  I.  Each  analysis,  involving  the  first  six  roots  and  associated  vectors, 
required  about  16  sec  of  CPU  time  using  the  UNTVAC 1 100/82  computer,  in  which  all  relevant 
matrices  pertaining  to  the  various  formulations  have  been  taken  into  consideration. 

Table  I.  Spinning  cantilever  bemm:  natural  frequencies  for  various  problem  types  for  a  spin 

rate  11 -0-1  Hz 


Structure  without 
damping 
(IPROB -2) 


Structure  with 
viscous  damping 
(IPROB  -4) 


Structure  with  viscous 
and  structural  damping 
(PROB-5) 


2-3953 

-0-3092  ±  2-3548 1* 

-0-3199±  2-35061* 

3-5689 

-0-3121  ±  3-54141* 

-0-3287±  3-53641* 

15-2142 

-0-3167  ±15-20771* 

-0-3912±15- 19641* 

21-7754 

-0-3166  ±21-77081* 

-0-4252±2 1-7643 1* 

43-0167 

-0-32021±43-01431* 

-0-5340±43-0298l* 

>  61-0008 

-0-32022±60- 99921* 

-0-6249±60- 98801* 

To  check  out  the  program  for  cases  VI  and  VII,  a  taut  string  vibration  problem3  waa  analysed 
using  the  EIGSOL  program.  The  results  were  in  very  good  agreement  with  that  presented  in 
Reference  3,  the  relative  solution  time  being  reduced  by  a  factor  of  about  2-3.  A  suitable 
problem  pertaining  to  the  quadratic  matrix  equations  (case  HI)  was  also  checked  out  by  the 
present  program.  A  cantilever  plate  free  vibration  problem4  of  aspect  ratio  1 : 2,  involving 
matrices  of  order  432,  was  also  solved  to  check  out  problems  defined  by  case  I.  The  solution  time 
for  the  first  six  roots  and  vectors  was  found  to  be  about  1  min  of  CPU  time,  compared  to  that  of 
2*3  min  using  the  program  EASI  of  Reference  4. 


CONCLUDING  REMARKS 

A  unified  numerical  procedure  has  been  presented  for  the  efficient  solution  of  free  vibration 
problems  of  usual  and  spinning  structures  with  or  without  various  forms  of  damping.  Such  a 
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formulation  also  includes  the  quadratic  matrix  eigenvalue  problem  associated  with  the  finite 
dynamic  element  discretization. 

The  program  fully  exploits  matrix  sparsity,  such  as  bandedness,  and  enables  computation  of 
only  a  few  desired  roots  and  vectors  without  having  to  compute  any  other.  In  general,  when  run 
on  the  same  computer,  the  present  program  is  found  to  be  over  two  times  faster  than  the  related 
programs  DAMP,3  EASI4  and  QMESSI.2  The  program  adopts  an  out-of-core  solution  strategy, 
and  as  such  it  is  capable  of  solving  large,  complex,  practical  problems.  Since  the  eigenproblem 
solution  time  is  proportional  to  NxMll3,  it  is  highly  desirable  to  adopt  a  suitable  bandwidth 
minimization  scheme  to  achieve  a  minimum  value  for  Ml  1  before  utilizing  the  EIGSOL  routine 
It  is  hoped  that  the  present  program  will  be  developed  further  into  a  «w»n  general-purpose  finite 
element  computer  procedure  in  the  near  future. 
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SOKMAXY 

This  paper  la  concerned  with  the  free  vibration  aaalpsla  of  coupled 
fluid-structure  ay* tees,  dlacretlsed  by  the  finite  element  method.  Both 
compreaelble  and  Incompressible  fluids  are  considered  In  the  analysis,  the 
litter  being  a  special  case  of  the  first  one.  A  numerical  analysis  proce¬ 
dure,  based  on  an  Inverse  Iteration  technique  In  eonjuetlon  with  a  special 
bisection  scheme  exploiting  the  Sturm  sequence  property.  Is  described  in 
this  paper  that  enables  computation  of  the  desired  roots  and  vectors  of 
the  vibrating  coupled  system  without  having  to  compute  any  other.  Further, 
the  procedure  utilises  the  associated  structural  stiffness  and  mass  natrl- 
e,i  no  well  as  the  fluids  counterpart  matrices  In  their  orglnal  banded 
form,  thereby  effecting  efficient  solution  of  the  eigenvalue  problem. 

XKTRODUCTION 

A  large  number  of  practical  structures  are  required  to  withstand  ex¬ 
ternally  applied  dynamic  loadings.  The  vital  preliminary  for  such  a  design 
requires  the  free  vibration  analysis  of  the  structures.  Involving,  by  far, 
the  major  amount  of  computation  time  of  the  entire  analysis  effort.  Many 
structures  exhibit  coupled  fluld-etructure  Interactions,  excellent  ac¬ 
counts  of  which  are  narrated  In  References  1  and  2.  The  first  category  of 
such  a  phenomenon  Is  characterized  by  large  fluid  notion,  an  Important 
example  being  the  flutter  of  aircraft  wings.  Some  numerical  solution  pro¬ 
cedures  of  such  problems  have  bean  presented  earlier3  **•. 

In  the  second  category,  the  fluid  Is  assumed  to  undergo  only  finite 
displacement,  the  notion  being  limited  to  small  amplitudes.  By  employing 
a  finite  element  discretization,  the  free  vibration  problem  of  a  fluid- 
structure  resonant  system  may  be  vrlten  as3 

K  m  -  m*M  u  -  C£  -  0  ...(1) 

H  £  -  m*0  £  “  *2P£r«  ■  0  ...(2) 

'where 

K  -  structural  stiffness  matrix 
M  -  structural  Inertia  matrix 

|  ■  fluid  matrix  associated  physically  with  the  Inertia  proper- 
~  ties  of  the  fluid  (analogous  to  R) 

0  -  fluid  matrix  asaodated  with  its- compressible  behavior  (anal- 
~  ogous  to  M) 

C  ■  fluid-structure  coupling  matrix,  sparse  and  rectangular  In 
“  nature 

u  ■  structural  nodal  deformation  vector 
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£  -  fluid  aodal  pressure  rector 
p  «  fluid  density 

end  in  which  the  fluid  idealization  is  aecoapllahed  by  the  standard  Euler- 
laa  pressure  formulation,  the  fluid  being  assuaed  to  be  invlscld  and  com¬ 
pressible  in  nature.  The  0  matrix  formulation  includes  the  coefficient 
denoting  the  speed  of  sous?  in  the  fluid  medium.  Equations  (1)  and  (2) 
may  be  combined,  as  below,  yielding  the  finite  element  fluid-structure 
eigenvalue  problem  of  the  entire  system3 


which  may  be  solved  by  standard  procedure  involving  real  matrices.  How¬ 
ever,  because  of  the  unaymmetrlc  nature  of  the  above  matrices,  the  solu¬ 
tion  tends  to  be  rather  inefficient  and  uneconomical  in  nature.  An  earli¬ 
er  effort6  succeeded  in  reformulating  equation  (3)  in  a  symmetric  form, 
although  the  entire  solution  process  still  required  a  considerable  amount 
of  computational  effort. 

In  Che  particular  case  when  the  fluid  is  assumed  as  incompressible, 
a  simplification  in  the  eigenproblem  formulation  may  be  achieved.  Thus, 
in  the  absense  of  the  0  matrix  an  expression  for  £  is  obtained  from  the 
reduced  equation  (3),  which  on  substitution  in  equation  (l)  yields  the 
corresponding  eigenvalue  problem 

[  |  -  ♦  pC  rlCT)]  a  ■  0  *••(*) 

which  has  the  effect  of  adding  an  additional  mass  matrix  to  the  structural 
eigenproblem  formulation.  The  final  mass  matrix,  however,  tends  to  possess 
a  rather  large  bandwidth  if  the  number  of  degrees  of  freedom  at  the  inter¬ 
face  happens  to  be  large,  which  in  turn  proves  to  be  expensive  for  the 
corresponding  natural  frequency  analysis.  A  solution  procedure  for  equa¬ 
tion  (4),  based  on  the  Inverse  iteration  method  is  presented  in  Reference 
7. 

The  purpose  of  this  paper  is  to  present  an  efficient  numerical  tech¬ 
nique  for  the  eigenproblem  solution  of  the  compressible  fluid-structure 
interaction  problem  defined  by  equation  (3).  The  procedure  starts  with 
the  natural  frequency  analysis  of  the  structure  in  the  absence  of  any 
fluid.  This  is  achieved  by  a  combined  Sturm  sequence  and  inverse  itera¬ 
tion  technique  that  computes  only  the  required  eigenvalues  and  vectors.  A 
special  inverse  iteration  scheme  ie  next  developed  for  the  coupled  system 
that  utilizes  the  eigenvalues,  computed  earlier,  as  starting  iteration 
values  for  convergence  to  the  required  roots  and  vectors.  The  solution 
process  cakes  full  advantage  of  the  relationship  in  the  relative  frequency 
values  of  the  structure  without  any  fluid,  structure  with  incompressible 
and  compressible  fluids  respectively.  America!  results  obtained  by 
solving  s  number  of  standard  test  cases  dearly  indicate  the  pattern  or 
root  convergence  corresponding  to  various  simplifying  assumptions,  furt 
demonstrating  the  relative  efficiency  of  the  present  procedure. 
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development  of  numerical  procedure 


To  achieve  as  effective  solution  of  Che  eigenvalue  problem  under  con¬ 
sideration,  equation  (3)  la  first  rearranged  as  follows: 


h  i  ol  { p]  To  IPC*1  M 

...»  -0 
-£  }  IJ  (uj  [o  I  kJ  [u  J 


which  say  also  be  written  as 


where 


(F  -  X E>X  -  0 


X  -  w2  ,  Z  “ 


•••(A) 

...C3) 

...(3a) 


Subsequent  analysis  procedure  Is  based  on  an  Implicit  assumption  that  che 

roots  computed  for  the  structure  (1^),  structure  with  Incompressible 

fluid  (X<X))  and  struct 
lowing  relationship9 •* 


fluid  (X<X))  and  structure  wit’  compressible  fluid  (X^)  bear  the  fol- 


»<»  >  >  x<c> 


..(*> 


'which  proves  to  bn  useful  In  the  determination  of  any  desired  roots  and 
vectors. 


The  entire  solution  process  consists  of  the  following  major  steps 
Step  I 

-  Solve  (K  -  XM)u  •  0,  the  eigenvalue  problem  pertaining  to  Che 

structure  only  to  yield  and  employing  a  combined 

Sturm  sequence  and  Inverse  iteration  procedure19. 


For  each  root  determined  In  step  1  solve  (F  -  »  0,  the 


eigenvalue  problem  for  Che  incompressible  fluid-structure  com¬ 
bination  by  setting  0*0  to  obtain  a  reduced  E  matrix  denoted 
by  K.  The  solution  li  achieved  by  an  Inverse  iteration  scheme 
and'a  bisection  strategy,  described  In  detail  later,  by  employ¬ 


ing  X<S)  end  as  the  starting  Iteration  root  and  vector  re¬ 
spectively,  yielding  sets  of  X^and  Xj 


Step  3 

Solve  (F  -  X^X>I)x  •  0,  pertaining  to  the  compressible  fluid- 

— -*  riX^  as  the  starting  Itera¬ 
te) 


structure  ease  by  employing  X^  end  Zi 
tlon  paraaetera,  as  la  step  2,  yielding  the  deelred  root  Xj' 
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and  vector 


Stop  4 


Triangular Ize  (F  -  '8)  and  perform  a  Stun  sequence  count  to 

determine  the  sequence  of  x'*C*  in  the  elgenspectrum  (say  rth). 

Then  perfon  a  bisection  procedure  based  on  the  strategy  as 
below:  tr\ 

(I)  If  r-1  then  X*  '  -  X*(**)  and  repeat  steps  2  and  3  for 
the  next  desired  root. 

(II)  If  r  >  1  then  repeat  the  Inverse  iteration  procedure 


with  X  -  (X<*>  +  X*fj)/2  or 


0S>  - 


2  (if  step  2 


has  been  omitted  from  consideration)  to  converge  to  the 

(C)  (C) 

•  required  1th  root  X^  and  the  corresponding  vector  ^  « 

Step  3 

Repeat  step  4  if  r  >  l+l  till  all  roots  up  to  the  1th  one  and 


Repeat  step  t  If  r  >  1+1  till  all  roots  up  to  the  1th  one  and 
the  corresponding  vectors  are  recovered.  .Assure.  AfCOcftW 


Step  6 


Repeat  steps  2  through  S  till  all  required  roots  and  vectors  are 
computed. 

The  Inverse  Iteration  scheme,  Implicit  in  the  above  steps,  is  carried 
out  by  utilizing  equation  (3).  Thus  the  iteration  at  the  rth  step  Is  per¬ 
formed  on  the  following  matrix  formulation 


r+1  „r+l_  r 

■  *g  lit 


which  may  also  be  written  as 


V‘tg | -v=Tl  fjTl .  -i  [ijLrjs'sL 

-£ ! i-'t"  J  [«ri  I  1  [  »«i 


...(7) 


...(8) 


where  lf^  is  a  normalising  factor.  Solution  of  aquation  (8)  is  achieved 

by  first  writing  the  matrix  equation  corresponding  to  the  lower  and  upper 
half,  respectively,  as  below 


(X  -  X^uJ*1  -  C  -  K^N  uj 
(H  -  -  xJpCTuJ+l  -  (S  zl  +  PCTuJ)nf'1 


...(f) 


...(10) 


The  procedure  starts  by  solving  equation  (9)  with  the  right  hand  vector  £ 
assumed  to  be  consisting  entirely  of  unit  scalars.  Equation  (10)  is  th*a 
solved  for  the  £  vector  and  the  process  being  repeated  till  adequate  con¬ 
vergence  Is  assured.  The  associated  root  Is  then  simply  computed  from  the 
Rayleigh  quotient: 


w 


vCv>- 
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I 


i 


»i  -  (zT'fciT1/  (zT1)1! iT1  — «*> 

Trlangularlzstlon  of  Che  relevant  matrix  pertaining  to  etep  2  le  ob¬ 
tained  by  setting  Q  »  £  in  equation  (8).  Indeed,  etep  2  nay  be  entirely 
omitted  for  compressible  fluid,  thereby  effecting  further  saving  in  solu¬ 
tion  time. 

numerical  RESULTS 

A  computer  program  has  been  developed  for  the  natural  frequency 
analysis  of  fluid-structure  systems,  which  is  capable  of  computing  the 
desired  roots  and  vectors  of  the  coupled  system  in  an  efficient  manner. 

The  example  problem  of  dry  dock  in  Reference  S  was  solved  by  the  program 
es  a  test  case,  and  the  results  correlated  rather  well.  Such  solution 
results  were  printed  out  at  various  analysis  steps  to  verify  the  efficacy 
of  the  bisection  procedure  adopted  for  the  present  analysis,  and  such 
results  confirm  the  reliable  nature  of  the  numerical  algorithm. 

CONCLUDING  REMARKS 

A  numerical  procedure  has  been  presented  that  proves  to  be  efficient 
for  the  free  vibration  analysis  of  fluid -structure  coupled  systems.  The 
fluid  le  assumed  to  be  compressible  in  nature  and  the  incompressible 
problem  is  only  a  special  case  of  the  generalized  algorithm  developed  in 
the  paper. 

From  the  numerical  formulation  depicted  by  equations  7-11  it  is  ap¬ 
parent  that  the  current  procedure  employs  the  individual  K,  M,  R,  0,  and 
C  matrices  in  their  original  banded  form  and  thus  fully  exploits  the  in¬ 
herent  matrix  sparsity  usually  associated  with  a  finite  element  formula¬ 
tion.  The  usual  procedures,  such  as  proposed  in  references  6  and  9,  in¬ 
volve  various  matrix  inversions  that  finally  requires  solution  of  eigen¬ 
value  problems  characterized  by  full  matrices.  A  similar  situation  is 
also  encountered  with  the  approximate  formulation  for  incompressible  fluid. 
The  p.-esent  procedure  also  enables  computation  of  a  few  roots  and  vectors 
only  without  having  to  compute  say  other. 
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